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Chaptei 1 
Introduction 


1 1 Review 

The held of tidnsmission and emission tomography has grown extensively m last two 
decides Dr Ilounsfield developed woild s first practical CT scanner in 1970 I his brought 
up 1 new u volution m (he fu Id of mi die d dngnosis xnd the techniquis lil« C ornpuUri/cd 
iomography (Cl), Magnetic Resonance Imaging (Mill) , Positron Emission Tomography 
(PEi), Single Photon Emission Iomography (SPECT), etc, are most widely used clinical 
and experimental tools Ihe success of these techniques lies in their ability to map some 
property of internal cioss section of the object without macroscopically damaging the ob 
ject ihis obviously leads to many mterstmg applications in medical diagnosis There are 
many other divergent fields such as radio astronomy , electron microscopy [1],[2] etc, which 
employ essentially the same techniques A common feature of all these problems is that 
1 set ol d xta is available that is related to the cross sectional properties m a known way, 
but not related straightforwardly The objective then is, using this set of collected data, 
to (stirmto tlu mxppiiig ol the piojicity iloiig the cross section 

In radio astronomy, the objective is to estimate the electron charge density distribution 
iiid migiKtu Ik Id iiil(iiHiiy didnhutiuii iiouiid IIk ou(<i stimiUK ol tlu sun [ Ij I lu 
(1 it 1 (olhctcd using tlu i ulio tc Us(op( iii ly is iii ( lu loini of s( i ip iiitc gi ils iiul tlu JD 
image leconstiuction problem is to be solved [3] 

Ihe tomographic techniques can broadly be classified in two categories, namely, trans 
mission tornogiaphy and emission tomography [1] In tnnsraission tomography the physi 



< il ()K)|>( I (y c V (linic <1 IS X i >y ibsoi pi ion ( 0 ( flu ic iil wliilc in ( mission I onio^i 
pliysu il piopcity cviluUcd is iiholon (mission inlciisily, ilso c ilkcl is iclivilion function 
Alllioii/ li ll)( l( ( lini(|ii( H of (hli (olidlioii iis( (I III 1Iks( (wo hkIIiocIh kc ( otisuk nhh 
(lill( i< 111 , I li( iiinj^i K ( oust 1 u( I ion piol)|( III I ill iiiidi I I Ik h iiik in il Ik m d k il loi iinlisin 

In CM scanning X rays arc pissed llirongli the surf ice under mveshg itiori to ohi un 
one dimensional projection at vinous angles which ire then sampled at uniform inter 
val llie data collected is normally contaminated with noise, owing to statistical nature 
of photon emission [2] 1 hese piojec lions ire eonvolvcd with a rirnp filler and the well 

known backprojection lechnuiue is applied on the filtered projections to reconstruct the 
(loss section 1 he i imp filtci knds to uiiplify the high frecpicncy components of data as 
well as noise, hut (he nois( usually doinuii((s il high litejucncics lienee (he r imp filter is 
typically combined with a low pass filter to avoid amplification of the high frequency noise 
11 k low piss fiKermg is also essentni (o avoid the iliasirig which otherwise lesults due to 
the finite value of sampling rate fioo much low pass filtering will result in loss of spatial 
resolution J he filters like Ram I nk. She pp J ogan,Gene rah/cd Hamming etc [4] are the 
most commonly used filteis for these applications 

Ree ently a ne w kind of low psss fillc r b ise d on llaar function w is dcsciibcd in [8] 1 his 
filter is angle dependent Another emerging mathematical theory which has found inter 
esling applications iii many engineering branches and in physics, is the theory of wavelets 
Wavelets have been applied successfully m various signal processing problems such as 
spcceh and image data compression noise lemoval edge detection ,ctc In tomographic 
ipjilu at ions wavelets have been used to construct small regions with reduced projection 
ehti [h], spue vnnnl noise filteiiiig [9 10] etc 1 he filter, when combined with the 
wavelet truisfoim based nonlinear space variant noise filtering have shown to yield better 
qu ihty of Ihe re constructed image This claim is verified through computer simulation 

1 2 Organization of Thesis 

• ( hiplei 2 d(s( iibes the basic piineiplcs of (i insinissioii loinogr iphy It ilso deseiibes 
in biief the two mam types of reconstruction techniques, i e transform domain and 
algebiaic technicpics 
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• Cltiptoi ] (Itsciihcs the lh<()iy lx luticl llu lint liltoi Tlu nlrvini pfitnlions irr 
tlso (1( 1 1 V( (1 III I Ins ( li iptt 1 

• ( *lnpl( 1 4 ^iv( s t 1)1 K I u ( ouiil ol 1 1)( \v tv( 1( ( 1 1 iiislotin 1 Ik oi y 

• ( 'lnpt( I 5 (t< s( 1 ilx I ( Ik ipplx il loll t (iI w \ v( I< ( 1 1 iinloi iii iii I oinoj^i ipliy 11 k ( iti 
pli ISIS will be on spUnlly vnyiiif, iioiiliiK u filUniig using w we let traiisfoiin 

• Chapter 6 presents the resuHs of computer simulation 
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Chapter 2 
Basic concepts 


2 1 Introduction 

Considei UiG space limited fund ion u{x,y) as shown m fig 2 1 We define a new 
fiiTution g{s 0) such lint 


/ oo too . 

/ u{x,y)6{a' cosO + ysm 9 — s)dxdy (2 1) 

-OO J — OO 

wltcre the ddt i function 6{x) h is following properties 

5{a;) = 0 , 5C ^ 0 



(2 2 ) 


/ OO 

6{r ~T)f{r)di = /(r) 

-OO 

(/(s,f;) IS ( died tlu H ulon ii insloim of i/(i ij) uul is dso Kfued to is proicdioii ol 
u{x,y) it an angle 0 J bus fi om eq 2 1 and the properties of S{x) we note that only 
Ihosc points (3,v) contiibutc to the integral which lay on the line P From hg 2 1, 
using the lotaUd co oiditnlc syslcin wc gel 


T = scosi? — 7/Siri d 


y s sm 0 u cos 0 


(2d) 
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i iguie 2 1 Radon transform as a line integral 

Using the above equation and the sampling property of the delta function (eq 2 3) 
wc get 


/ “ 

= / u(a(os<? — u sin 0, s sin 0 + ucos 6) du (2 4) 

J — CX) 

In practice 1 heso ray iritcgials aie obtained by passing X rays through the object un 
del tesl md Ihe lunctioii u{x^y) rcprescnls the distribution of the X ray absorption 
coefficient (In case of PUd k SPPCr, u{x,y) represents the activity distribution 
[1] ) Iho pioblem is then that of inversion or finding the function u(x,y) given the 
piojpctions g{'i,0) for all values of 0 in the range [0 tt] 

2 2 Discretization of the reconstruction problem 

In practice the projections are obtained over a discrete set of angles and each projec 
tion IS sampled at equally spaced points Let gn{m) denote the sampled projections 
Jhc discrete piojcctions arc thus defined as 
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m e {-M -M + l -L 0 1 M] 

Hul IS I he s inipliiig I c III pio|( ( I loii <loin mi mcl is sc Ic c t c d such { h \ l 



(2 b) 


D specifies Uic boundary of the region wiUnn which g[s 0) is defined i c (j(s d) = 0 
V(? if |s| > D Also the estirm.tp of (he function v{i y) is e\ dilated at disciclc 
vihies of 7 &; 1 / which iic spared it unifoim interval Two nielhods aic used foi flu 
computation of u{x,y) algebraic reconstruction techniques and tiarifaform doiiiiin 
nictliods 


2 3 Algebraic techniques 

let us assume lint (ho function 7c(r ly) can be oxpie sse d by i finite series of the loiin 

p-i 

y) (2 ' 7 ) 

=0 

whore aio known functions, typically reetangiihr functions with noiiovcil ip 

ping region of support as described in [6] The coefficients {y } in the senes expansion 
,the n elnt lete 1 1 /( (lie tune t luii i/( /, ly) e oinple (c ly let /fj denote (he Rulon Inns 
form operator conespondmg to the piojection simple e/j Using the Inuanty 
jiropeity ol the i uloii (i nistoim we get 

rtAi,y)=^ Y.v RjJ [x y) (2 S) 
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'/ j = Z! !/ “ j 


(2 9) 

where 

« j = (2 10 ) 

Let G denote the eolumri vector of the sampled projections and let V denote the 
vector of coefficients in the senes expansion Then the problem can be modeled as 
that of solving the system of simultaneous equation 


G^AY (2 11 ) 

The mitiix A IS tjpually a lectanguhr rnatiix with vtiy large size ( Of the order 
of 10 YIO’’ or moie) This makes it virtually impossible to invert the matrix A and 
lieiKC itciative techniques ate used lor the solution of the vector Y The vector Y 
IS given 1 lint 111 vahu and at ( uh iteration some correction term is applied to Y 
so tint the vedoi Y convergis to, oi is dose to tlu solution of equation discussed 
m [6] iheie are various types of iterative techniques such as additive ART and 
miiltiplK ilive ARi These imthods differ m the manner in which the corrective 
tcirn IS qiiiludtoV [G] 1 he m iiii dr iwb uk of the iterative techniques is that they 
deni iiid much c ornputational load and the rate of convergence is very slow Hence the 
Algebraic reconstruction methods are rarely used although they give better quality 
of the reconstructed image than the transform domain method when the no of 
projection angles are few or when the projections are incomplete [6] These kind of 
situations oiciii, for exmiple, m fist cardiic irn iging or in 3D image reconstruction 


2 4 Transform domain methods 

All of the Transform Domain methods are based on the inversion formulae developed 
by John Radon m 1917 The following conventions are used in the discussion to 
follow and in siibscquciil (liqilcis 

Small letters (eg g,u etc ) are used to denote the luiictions (of one variable or two 
vunbhs) III spue doiiniii iiul c ipit d Icllcis iic used to dc note tlu fmictioiis in 
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licqucncy domain llius G{(,0) dcnolc one dnrieutionil founer transfoim of g(s 0) 


/ OO 

g(s ^ ds (2 12 ) 

-OO 

1 Ik two dinunsioinl louiKi (i uisfoitn of v(^i ?/) is similarly defined as 

/ CO rca 

/ v(t + dxdv (2 13) 

-ol J xj 

uid the irivc rs( loiiru r transforms nc given by 

/ OO 

(2 14 ) 

-OO 

and 

/ OO roc 

/ //(6 (4 15 ) 

-OO OO 


Wc denote the polai representation of a function of two variables by putting a hat 
on it Ihus 


u(r, i/;) = u(r cos 0, r sin i/j) (2 16) 

f/(^,(^) = (/(^ cos (/!>,{ sme^i.) (217) 

J hiis the Fouriei transform and the inverse Fourier transforms can also be expressed 

IS 


roo rlTV 

U((J) = J j (218) 

u(r,i,) = /“ J‘' (2 19) 

For the function f/(iC, <?i) we have (lie symmetry relation 

II </,_7r) -=/;(( e/i) (2 2Uj 
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Using tins icKlion it cm be shown 111 it 


/ CO /•TT 

/ (2 21 ) 

•<x> 10 

Accoiding to the ceiitial slice tluoicni [4] wc get 


({t, 0 ) - {J{iO) 


— //(iCcosd ynid ) (2 JJ ) 

combining the h suits we get 

/•TT roo 

u ( r , y )= / / ( 223 ) 

Jo J- JO 

lutlu xbovc cc[u itiou the iiinc nnost intc gt il w i t i( can be intc rpreted as hltcrmg the 
piojections with a ramp filter while the outer integral w r t t? is the backprojection 
opciation Let /3 denote the backprojection operator Thus we get 

ii(] — f (j(t 0 + y sw 0 ) dO 
Jo 

ihe ramp filter tends to amplify high frequency noise and is unstable, hence it is 
typically combined with a low pass window function VU(if) and the inverse relation 
be comes 


/*7r /* oo 

u(x,y) = / / 144(0 (2 25 ) 

Jo J-co 


= l^Pg 


(2 26 ) 


whcic 


/ OO 

r ( 00 f c/e (2 27 ) 

-oo 


= q{SjO)*w{’i) 
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Figure 2 1 Ireriuency response of different convolvers 


iu[b) IS tlu OIK dime iiHioii il iiivtise fouiui Iruisloim of 

Various windows have been proposed m literature and the most widely used windows are 
Ram Lak, Shep I ogan, Generalized Hamming and Low pass cosine The plot of 
tor ill these window functions are shown in fig 2 2 


2 5 Discrete Implementation of Inverse Radon Trans- 
form 

The impulse response of the convoher w{t) is obtained by inverse transforming and 

tire discrete convolver is obtained by sampling the continuous version with proper sampling 
rate [4] 

Let w^w) denote the discrete convolver and let us denote the filtered projections by ^n(^) 

Then ^ 

9nim) = ^ gn{m - k)io{k) (2 29) 

k=-M 

As the function is defined for integer values of in, the value of the projection at other 

values of m is evsluated by iiiterpol ition I inear interpolation is widely used in practice [4] 

1 bus ki p, d( note the integer pari of a red number p ind let pj denote the fractional part 
of p Ik nee p E Z and 0 < 71/ < 1 1 01 linear mterpohtion ,the interpohtcd projection 

9n{p) derined as 


11 






Pol 11 C 1 1(1 



He ( laiigiil u ( ncl 
1 igut( 2 3 C nd SLiu( tuic s 


= ( 1 - {i)i ) \ rn j (j [] -{■ ni ) {> JO) 

1 Ik 1) K k|)io)(( (ion opi i iloi is ( lu n disculi/cd uicl llu disc u ti/(d f slim ilc of (lu [niulioii 

i/( I , I/) IS ohl (iiiK d IS 

N-l 

/ f I \ TT . yllTC 

“ 77 E •Mf' )y ) I /«i'i(- )) {> ji) 

wluu ~ IS ( he s impling i lU iii sp ill il doin iin uid (L 1) £ 

llu Him 1 ak (iH( 1 II IS sli II p ( 111 oil Ihis usiills in high o\( isliool uid i iiiging in I In 
region ol high contiasi i he ollui fiKds h iv( coin]) iral iv( ly snioolh loll ofl 

2 5 1 Fourier Domain Reconstruction Methods 

AnoLhci li nisfoini doinun Icchnique no! so widely used is the fouiui domain rcconsliiu 
tion method [7] In (Ins method the projections aie fouru r liansfoi me d Ihisgiv(s (lu 2D 
loniu 1 (lansloiin ol i/(i , i/) sainiih d on a polai giul is shown in hg 2 J I lu pol u gi ul d il i 
IS ( on\( I (( d (o 1 i( ( ( ingul u gi id d il i using iiitc i pol ilion ind ( lu 21) iiive rsc 11 1 is ( ike n 
so IS (o ob( nil s irnplcs ol 1/(7 7 /) over 1 ier( ingul ii giid Vaiious in(c 1 poKtion 1( e hiiujiu s 
ue piojiose d (0 (onv( i( pol 11 giid (hi i lo iee( ingul ir giiel [7] llu m(e tpoh( ion e iioi is 
(he nnin souicc of eiioi lienee hltcicd baekprojcction tee hnic|iic is jiie fi tie d ove i fouiui 
domain leconsl ruction methods 
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Figure 2 4 Parallel beam and fan beam geometry 



Figure 2 5 Fan beam projections 


2 6 Fan Beam Reconstruction 


The reconstruction geometry , considered so far, is called parallel beam geometry Another 
me t hod widely use d in practice due to it’s siinphcity is the fan beam geometry The elet uls 
are shown in fig 2 4 

In fan beam geometry, the source moves m a circle of radius R, while the object is m 
concentiic circle of radius D/2 The detector is either in the same circle as that of the 
source oi m a circle of diffeient ladius, and is diagonally opposite to the source as shown 
m fig 2 5 

The cone beam projection is denoted by b[a^P) where a € [—u u], Q < I3 < Itt and 
V = aicsin(^) 
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A fan bcTTn pro|fclion dxta is ilways some parallel be nrn data Thus fiom fig 2 5 we got 


/j(cr, j3) = g(R sm <7, er + /3) 


(2 32) 


9{s 


6) = 6(arcsm — , 0 


arc sin — ) 


(2 33) 


Thus the parallel beam reconstruction method can be applied This is called as method of 
rebivnirig Diicct reconstruction methods are also possible [3] 

A squaie image and the Shop 1 ogan Phantom [11] are used in simulation Ihe results 
xre picscnled in chap 6 
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Chapter 3 
Haar Filter 


3 1 Introduction to Haar filter 

Ihc fillers rnf niioucd in cli iptf r 2 irc all band liiniled (m other words they are all having 
finite region of suppoit in ficquoiuy domun) ind hence they can not be space limited 
The band limited filters allow one to avoid the aliasing when the sampling rate is properly 
choosen so as to satisfy Shanon’s crilcria, bul in praclicc the projections are not low pass 
filteied prior to sampling (although the sensing device provides some sort of low pass 
fill (ring, it IS not adecjualc), hence aliasing is always present to some extent in the sampled 
projections On the other hand the Ilaai filter has finite region of support m spatial 
domain resulting in some aliasing, which can be contiolled by oversampling ratio p (see 
next section) Ihe overall effect is a belter quality of reconstructed image 
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0 * v ) 


y 

1 i{,ur( 5 1 21 ) II iir (iindion tiicI it,’s pu |(clioii 
(\)nsKl(i the 2D 11 ur function shown in hg 3 1 1 hifa fuiuiion is clehncd as 


•So A(a",y) = Soa(-c)so A(f/) 


(3 1 ) 


111(1 flu (uiulion S()a(’) is(I(Iin((l is 


So a(2) 


1 il \z\ < Y 


z = 


{ 0 if 1^1 > 


I ol u(t,7/) be a continuous function in L"^(E) aiui let Vq clcnote the space of functions 
which IS spaiicd by puccwise continuous luiictions Sq ^ — dil) wlicie {LJ) € Z'^ 

1 ct u{x,y) be the pro)cctiori of u[x^y) on Vo Ihen we get 

CO oo 

n{i,y)= X) Wo(A-,/)so aI^- - AA .,1 - A/) (12) 

fc=— oo /=— oo 


I Ik f-i iiiiph s u( givinby 

Uo{k, 1) = ((u(a;, v)), (^o A(a: - ^k,y - A/))) 


(3 3) 


where 

((^^(^, 2/)), (so aI-t- — AA,, 7/ — A/))) = / / u(a,,T/)Ao A(a — AA,,i/ — A/)dr(/y (3 1) 

J—oo »/—oo 


1 ') 



|(l /^(Utu)(< (lu t uloii I ! inslot in ti|)( mIoi nnl I< I // (It ho(( IIh !> n kpiojt ( I hjii opt i ilor 
IS disc usscd m tlip 2 1 loiii Uu i)i( vious discussion the lollowiii^ relations iiold 


UQ[k 1) =^u{x,y)\^=cky=M (3 5) 

^(6,6) = f/(6,6)5oA(6,6) (3 6) 

Let l\i denote one dimensional louiier trmsfoim ope r dor md F^d denote the nivers 
feniiiei 1 1 nisfenin 0 [)e i dot We define some piojeeiions itid ilitii It insfoirns is in e e| 17 


to eej iliesc lelations aie used in further diseussion 

g{s,9) — Ru{x,y) (3 7) 

g{s,0) - Ru{x,y) (3 8) 

uJoa(‘ 5 (?) = iiso A(ai,y) (39) 

Wo^[i,e) = rriWo^{s,0) (3 10) 

C[U) = Fid9M (3 11) 

G{(,0) = ru9{s,0) (3 12) 

Using eq 3 6 to eq 3 12 we get 

= (3 13 ) 


Ua{k,l) = /l(|^lfUo(^j ^)^u(i^) <^))|r=AA: v=Al 

= l3K{^,e)Gu{^,0) (3 14) 
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where 


= 0 ) (? r,) 

Thus the filter k{s,9) = becomes the reconstruction filter in the convolution 

bdckprojection algorithm The sampling kernel So^{x^y) along with its projection is 
shown in fig 3 i As the sampling kernel So^(a; y) is separable in spatial domain, it is also 
separable m frequency domain, le 


VU„a(^,^) = H/o(^ sill C>)VKo(^ cost)) 

{S 16) 

l'ko(<C) = 1 lrf‘>0 a(^) 

(3 17) 


U can be sliown that 


Vko(0 = Asmc(eA) 


(3 18) 


Hence we get 


IKo a(^i " A'* Mine (A^ (OS 0) siik (A sin 0) 


(n<)) 


A (i^, 0) = |<^|A‘* siiic(A^ cos 0) sinc(Aif sm^) (3 20) 

Ihc (xpiessious foi k{s,0) = if/A(<^,^) arc dciiv(d in appendix A Using these results 
we gel 


fc(s, 6 ) 


(b - ^(1 + s,n 2«)| - In 1 5= - ^(1 - 3m2«)l) 

, 0^0 


k{s,0) 



This IS the expression foi the reconstruci ion filler m the continuous case 


(3 21) 
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1 iguic i 2 Plot of VPo a(^ 0) 

3 2 Discrete Haar filter 


ihe discrete haai filter kiirn^O) is obtained by sampling the continuous version Looking 
at the plots of Wqa{^, 0) (see fig J 2) it is clear that if the sampling interval r in projection 
domain is laigcr or even equal to A, strong aliasing will occur for all values of 9 The 
resulting artifacts degrade the reconstructed image consideiably Hence the projections 
are oversampled by a factor p > 1 and 


A 

r 


(3 22) 


Also the function k{s,0) described by eq 3 21 has some singularities at (a = 0,0 = vr/d) 
and (s = A/4,0 = 0 5 arcsiii(3/4)) These singularities occur at sampled values of the 
trapezoidal shape where the derivatives are not defined Hence the filter k{s, 0) is integrated 
by rectangular window and then sampled to obtain the discrete filter Thus the projection 

IS sampled at 


and 


s 



P 


(3 23) 
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k i{i} 0) = 0)*i((l{ — )^ 

/■{ + i ) 

= k[uO)du (J21) 

1 lu i)io)( ( 1 ions i()|) a( ■’ 0) IK ( (I lo 1)( s[)( ( ifu (1 foi 0 C [0, 7r/-1] 1 oi ol lx i viiliu h ol 1) I lx 
folUnviiif^ lol ilions hold 


xj,) a(s 0) = lOo A(s,7r/2 -(?),(? e [tt/ 1 7r/2] ( 5 2h) 

<"()a(s d) - t»oA(s — 7r/2) d 6 [7r/2 Jtt/'I] ( 5 

"’ll a(s d) - ID,) a('’ TT - d) d G [ Stt/I H-] ( 5^'?’) 

ll(ix( llx luixlioii Iho A((,d) lollows till siiiK syiniiuliy locviliiili llii (IihikU fillii 
W( iioU 111 il 




ds 


[(s - r) 111 |s - r I + (r + s) In ji + s|] — 2 


(i2S) 


'Thcicforo integrating the above function wc get 


kAn 0) = — [(/ — o) In |/ — al I- (t + fv) In |/ + rv 

TT sin 2w 


-(/ - /^) In \L ~ ft\ - (/ + /i) In |/ + fk\]\ OjkO ( ! 2<J) 


W lu K 


a = — a/ 1 + faiii2t> = ^(cos + sill 0) 
2 2 


(i 50) 


§ — — \/l - sill 2Q — ^(cos 0 - sin 0) 
2 2 


(i 51) 


23 




1 igim \ 3 fcrf(n 0) p ~ 2 

I he ibovo expression (or tj(n, 0) ,s vn[,d for 0 / 0 h{n, 0) ,s similarly obtained by filtering 
and discretization of k{s,0) Thus 

Arf(n,0) = — ^3 32^ 

lli( list ((ludion cl( lily shows lint tlu fdtd ki{n,0) is still not d( fmc d foi odd vilucs 
cl /I which Ins singuhi.tus it n = ±(1 + /i)/2 uid n = ±(1 - p)/2 1 his is the major 
limit ition of Ilaxr filter 

1 Ik [dots of Arf(7),d) for vinous viluos of 0 ^Ild p nc shown in fig 3 3 to fig 3 10 
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0 1 1 \? 



n -4 

i igun A:/(n,7r/]2) p = 2 


9 — 7r/6 



1 iguic 3 5 kd{n,'K/C)) p ~ 2 
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n 


Figure 3 10 A:d(n, vr/d) /? = 4 
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Chapter 4 

Wavelet Transform 


4 1 Introduction 

loi loiK (Imuii loinl lun(lK)n7(,) Hu lontid It uislotrn / is uioihc r k prc sc nUl ion of 
llu sunc liiiu t ion s|)( ( ific (1 111 Itcciitc luy (loin uii 1 lie louiu r 1 1 iiisform (duel subscc|ucatly 
tile inveise fouiicr iransfoim) is oithogonal and all information m the original signal is 
pieseived in the liequency domain Note that we can write the inverse expression in the 
form 

/ CO 

< dx (4 1) 

Ibis indicitcs tint the function f{x) is expressed in terras of it’s projections on the basis 
functions Flu sc functions aie peifectly localized in frequency domain, but have in 

finite extent in tune doin nn lienee iiiy time locdi/ed inlornidtion in f{x) such is ibiupt 
changes ,is spu ad in the whole spectrum Hence the fouricr transform larks the knowledge 
sbout the lime loe ili/ed (catuies ol the sigiul lo address this problem, Caboui mtro 
duccel the window fomici transfoiin, also called as ’short time fourier transform’ (STFT), 
or Gaboui trinsJoim in 1945 

The S IFI of a siguxl /(/) is defined as 

/ OO 

(4 2) 

-OO 

where li[t) is s low psss window like Canssion window, which is locah/ed in frequency 
doin 111 ) ns we II is in I line dom iin 1 lie {> n mu te i (, spe e die s the eoiive nt ion il fie (pK iic y incl 
T specifies the time shift Thus the ti insfoim SJ TT;(Gr) has significant magnitude only 
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]f tlip signal /(f) his significant spcctial component of frequency ^ in the neighborhood of 
llu (iiiK insl lilt / — T I( IS shown llnl logel fiiu Ik (|ih lu y u solution , the window should 
he IS wiele is [leissihle , iiielllie c eiiielil ion of (me lime le solul ion ele iti iiiels sin ilk i siqipoi 1 
veiiidow I li( i( ( oiiIIh ( iiif i< (|iiii( me iiIm < III iiol ohvioiiiily he ii il mile <1 iimiiill me oimiy mel 
llie IJnee il iiniy Piine iplc gives the lowe r hiriils on Ihe pioelue t of ihe time scale resoliilion 
and liequcncy scale resolution [18] ihe SIFT uses a window of fixed support, hence of 
fixed bandwidlh This type of analysis is called uniform filtei bank inaly&is [12] 


111 this regaid the Wivelel transform offers the advantage over the other time frequency 
elisl 1 ihiil lems It dleiws one to hive (me lime se ile Ksohiiiem it the eosi of fieepieiiey 
scale resolution and vice versa (Within the limits imposed by the Uncertainty Principle, 
e)f com sc ) Along with Wavelets eorries the notion of scale The continuous parameter 
Wavelet ti irisloim is defined as 


1 X — T 

Wj{t, a) = I /(s)0( ) dx 

■\J CL J —CO (I 


(4 3) 


ilie window function i/i(ec) , called as mother Wavelet, is a bandpass function In Wavelet 
tiansfejiin, the analyzing window is scaled (le either dilated or contracted) and the shifting 
parameter t is made a function of the scale parameter a Thus narrow windows axe shifted 
by small interval and wider windows shifted by large intervals This results in ’constant 
Q ’ hltei bank analysis [12] Ihe Wavelet Iranslorm given by eq 4 3 is invertible provided 
the mothei wwelet satisfies the admissibility condition [18] given as 

l^(OP 


/: 


iC 


• de = 1 


(4 4) 
(4 5) 


whcie 

\y(i^) = / dx 

J —oo 

When il){x) is having fast decay at mfiiufy [18], tin ulmissibihty condition can ilso be 
wiittcn as 

'ip{x) dx = 0 (4 6) 


/: 

Ibis indicates that tjj{x) is oscillatoiy 


4 2 Discrete Parameter Wavelet Transform 

The Wavelet transform given by eq 4 3 is very much redundant The same is true for 
the STI d also Hence it is desirable to discretize the parameters (a, r) ((if, r) m case of 


10 



T 


I ijjUU 4 1 Dyitlu ^rid lot disc ic 1 1 / i( ion ol («,r) 

SIl r ) In fact we would like to have a specific grid for (a, r) such that the associated 
wavelets {ij) t( 0) orthogonal wheu 

= {- 1 7 ) 

y/a a 

1 he orthogonality implies no redundancy It is shown that for SI FT ,wh.en the sampling 
gild is selected so as to aehieve the orthogonality, the corresponding window functions are 
vciy b idly localized either in time domain or in frequency domain [15] This is why most of 
the tunc tlm S J F 1 is oversamplcd IIowc ver foi Wavelet transform , it is possible to achieve 
oithogon ilily along with loeali/ation [I*!] 1 he discrete parameter Wsvclet transform is 

defined as 


DWTj{m,n) = lF/(a^,n6oa^) 

- n fOO 

— / /(x)7/)(ao”'a: — nfco) dx (4 8) 

J —00 

Mayer has show that there exist a class of wavelets {•0 (t)} for which the choice cq = 2 and 
ho = 1 gives oithogoiid hiscs of /^(/O ••'is concspoiids to the dyadic gild as shown in 

hg 4 i 

4 3 Wavelets and Mult iresolut ion Analysis 

The multiresoluiion analysis methods are widely used in computer vision for the purpose 
of pattern recognition In this section we will consider in brief the relation between 


II 



mull iH soliil KJii 111 ilysis iiul w iv( Ic ( liuisloim A iiiok (1(Iu 1((1 lie iliruiil c m Ijc found 
in [12] 

llu inuKiicsolutiou analysis consists in bi caking up into a ladder of spaces Vj 


CViCViCVuC K-i C l/_i ('i ^1) 

11k oithogonal piojc c I ion of soiTK fundion f{r) £ on V, rorrosponds to approxirn i 

tiori ot J{'i) at icsolution 2^ the Mnllnosoliilion ippioxunation his following jiropcrtics 

- fi2^x) e Vo 

“ There exists a scaling function (^(-c) such that {<po n{^)}nez is orthogonal basis of 
Vo wheie (j)on(^) = <^{x—n) This also implies that = 2^(f){2~^x — n)| 

,n e Z forms the bases of Vj Note that the factor 2^ is present so as to make 
the norms of (ji{r) and i(r) identical 

- Let us defiiK the orthogonal compliment of in by Wj ,ie l/j_i = 0 kK, 

Ihen theie exists a function i/>(x) such that {i/’jnl'?')} 's orthogonal bases of W^j 
where 

il)j^{x) -2~^i){2~^x - n) (4 10 ) 

— I ct Aj denote the orthogonal projection operator on We denote the projec 
tion of J{x) on V, by AJ{x) Then AJ{x) is completely characterized by the 
sequence {n^n} where 

Cijn = < > (4 11) 

- I he approximation at a resolution level j is coarse approximation as compared 
to approximation at resolution level j - 1 d he resulting loss of information is 
given by orthogonal projection of f{x) on Wj and is called as detail signal at 
a resolution level ] Let D, denote the orthogonal projection operator on the 
subspace Wj Then DJ{r) is completely characterized by the sequence {d^n] 

where 

djn = < /("c) '^jn(T:) > (4 12) 

— llie functions <j}j (r) & Vh«(^) satisfy the following relations 

< <t>]k{x) , > — ^k-m 

\2 



J-l 


3— m£)~* k2) 


L- l » I — — ^g)- 


1 1 

Uh^ 

1 1 




C 1 

uJ 



j-i 


I i^UK 4 2 idiuuul piclcd k (oiisl i lu tion (]u idi dun miiroi fillci bank 


< Vbit(3) <f>i (t) > = (4 13) 

lioin the foicf^onig discussion d is clcai that {?/>j J ,(n,;) G is orthogonal basis of 
Mallat has shown [12] how to obtain the sequence {ct^ n] i^jn} given the 
sequence „} and vice veisa This scheme is shown in fig 4 2 

Ihcse sections are identical with the analysis and synthesis sections of 2 channel , perfect 
i( construction quadrature mirror filter bank (QMr)[19] The filters, G, &:// satisfy the 
conditions 


l/f(a;)p + |//(a; + 7r)p = l 

(4 14) 

II{z)Giz-^) + H{-z)G{-z-^) = 0 

(4 15) 

Il{z) := X: 

(4 16) 

TL^Z 


(4 17) 

nG7 

Cr{uj) = G(/)| . 

(4 18) 

H{U>) = 

(4 19) 


Ihe concept ol multuesolution resolulion can be extended to 2D images [13] 






I igut( 4 ! Pyiuiiulil ilgonllmi Ati il> sis section 



1 iguu 1 1 l’iiutii(Ul Ugoiitliiu Sj iitlu SIS sc ( lion 


4 4 Fast Wavelet Transform 


In |)i u 1 1 ( ( I Ik sigrnl f{i) is s unpic cl willi firiik s iiiipling i ilc I / / , to obi uii I In sc epic iic c 
{>i{ii)} 1 he sc s iiiiplc s c in lx suppose cl lu give pi o)c c I ion ol /(f)oiil()i( 

(yQn=u{n) I he wavc’lci c ocflicic Ills c/j I > j > 1 ran be obi aiiircl using the pviiiiiichl 
nlgoiithni IS shown in (ig 1 i 1 be synlhcsis sclieiiie is shown m hg ^ 'I Ibis is llie I isl 
wive let It uisfouu ilgoiilhm 


4 5 Compactly supported orthonormal basis 

Using (he eonapis dciclopecl by M ilhl Diubeehies lias shown how lo obi un eoinpielh 
suppoited wave 1(1 bases bom Ihe (|uuUduic minoi hlleis [20] J he Diubeehies wive lets 
have finite tegion of suppoi t 1 he length of llie fillet /i(ti) \ tiics aecoielingly I Ik plols ol 
wave lets iiid se ding luiiclioiib lot diflcie ill fillci Icnglhs aic shown in fig 1 5 lo fig 1 '3 J he 
filte 1 eoeffieie ids ite give ii in [10] I die is of ie iiglli d 8 iie e oiisiele le c| (oi ele nioiisl i ilioii 

llie oil hogoiid wave le Is mid se ding fuiiclions mgeneid ue not syniineliie Ihisiesiills 
in nonline II plnse fiKe i s in inil>sis iiicl synlhcsis filter Innlcs By ulopling bi oil hogon il 
fillets ,thc condilion ol line at phase t ui be salisficd However this is at the eost of using 
sepaiale sc ding fuiulion ind wavelet in inalysis and synthesis sections ol Ihe niultiiej^o 
lutiou dceoinpohilion 1 oi addiLioiial etc lads the re idei is icfciied to [17] 
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Chapter 5 

\Vavelets in Tomography 


5 1 Review 

As mcnlioncd in the hist chdpter, the lime frequency localization property of wavelet 
t lanslorni is ul ihzcd in I he applications like region of interest tomography [5], space variant 
riois( hltonng [9] ,[10] etc In icgion of intoic&t tomography, only a small cross sectional 
11 oa IS Kconstiuetid at a highei u solution Using wavelet expansion of the image ,th]s 
iccoiistiuction can be carried out using less no of samples per projection [5] , with a 
subsequent reduction m X ray dosage Another application of wavelet transform is the 
spall illy V 11 ying flit CIS for noise sujipKssion Once again the space frequency localization 
pioperty of w \v( 1(1 ti uislotni is nude use ol In tins ch ipter this ipplic itioii is considered 
m more det uls 


5 2 Space varying filtering 

J he piojections sampled in the tomography are normally corrupted with noise The stalls 
tical nature of photon emission is the mam source of noise [2] This noise can be modaled 
as additive white gaussian noise with zero mean [11] Thus we get 

9r(S^) = + //(s t?) (5 1) 

where <7, (s,d) is the corrupted projection data, 7 (s,t?) = Bu{x y) and w{s^0) is the 
additive white gaussian noise with zero mean and van ince The noise component is 
uncoi elated in both the variables s and 0 Thus using the expectation operator [22] we 
get 



I [/'(•’I 


Uhiii/ llu H uii|)l( (1 v( I sions w( {,( t 


(5 2) 


9n(w) = + /2„(m) (5 3) 

I oi llic sake ol sirnplirUy wo will drop Iho lerm n is the forogoing discussion is ocpiilly 
ipplit il)l( to ill (he piojcctioii dugk s 

] ci d j d j and d ^ denote tlic discrete wavelet Iransform coeflieierits of noisy projections 
, uncoil upt ( d piojcctions and white noise component icspectively, at resolution level i 
Using the line inly propeity of (he wavelet iransfoim we get 


For icduction of noise power ,thc pro)ections can be low pass filtered using linear filters 
Plus liowtvti , results in loss of spatial resolution The solution proposed to this problem 
[5] IS a space van int filtei in which, only those wavelet coefiicients, which have magnitude 
below a ccilsiii (hicshold, ire set to /eio As the wavelet coefficients represents the detail 
signal 01 high ficqiuncy components, this is equivalent to space variant low pass filtering 
1 lu piocesscd wwclcl coelhcients d,j uc thus given by 


d 


1 3 


■ d:, ,f Kj > A{.) 
. (I ■! < -'(•) 


where \(*) is some I hicshold corresponding to level i 


(5 5) 


5 3 Selecting the threshold 

Foi threshold selection we use the energy cniena given m [24], [23] 

The threshold A(?) is selected such thit the total erieigy of the processed detail signal at 
ic solid ion level i is close to the total cncigy ol the wavelet coefficients d ^ (which corre 
sponcls to imeoirnptcd pio|cetioiis) d (he same resolution level For orthogonal wavelet 

lepic sentalion we have 
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jV (5 6) 

3 

wi„K ,hi . ,i v„s.,i„i,„n 

|( V< I / 

I x|) uidiiifi, d)()V( < <111 iliou w( }i,( ( 


1 

(5 7) 


(5 8) 


J J 3 


Note llul lh( wluk noise is sssuTned io bo umonhUcl with g (s,^?) Tiid is also assumed 
lo be of /( lo me ui md van uicc <r^ lienee wo gel 


EK/,,] = 0,V^,VJ 

(5 9) 

n{<,r\ = 

(5 10) 


(5 11) 


lie re h ele Holes I ho expeeialion oporaioi and it can be shown that [22], [24] 


)(E ^ (5 12) 

n n 

Jleiue fiom eej 5 8 &■ 5 11 we get 


/ 4 /^;] = P“d h[P''] 


(5 13) 


wlu re 


h) 



MP'l = riU'f:?] 

3 

- LMW',)'! 

1 

3 

I(t M dcnolr the lol t 1 no of s miplos jx i piojt < ( loii As Uioic is i dc ( ini iLioii !)> iluloi 
of two xllci evGiy filtciiiig sLagp of the anily&is section (see Rg 4 3) it follows tint i( Ai is 
the no of wavelet coefficionls at usokition level i then Af = M iiul 

M -I 

'l'"l - 

3 

- I\l2~frf () 1 ^)) 

Liu Ihushold \(z) IS sc luted such 111 d ill 1 he v\ ivc let coe Ihe lent s ( il a lesolulion h \c 1 
/), whose inagiiitude is below the threshold hive i totil cneigy iiiinoxiin ite ly e<|Uil to 
I [Pi'] 1 lu pioc ( ssc (1 w ivc 1 ( I e o( file K Ills nc Iheiiinvcisc 1 1 ilisfoi me <1 iiidllu pioecsscd 

piojcelioiis ue used foi linage ic eoiisti ucLioii 

J Ills lcchnu|uc has aheicly been used siicrcssfully foi noise hltering uid is shown io give 
bettei pcilouiunce conipircd to linear space invanant filters [23 9] Llie objective of using 
this method hue is Io see the combined effect of Harr filler and wavelet pic piocessiiig 
1 he lesulls of eoinpulei siniulitioii iic give ii in ne xl ch ipte i 
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Chapter 6 

Simulation results 


6 1 Staiidard Reconstructions 


In tins (liLpiti wc piesent the results oi computer simulation to demonstrate the elFec 
tivciKss ol the IHu filter In the first section wc consider the reconstruction of a square 
inn^( I h( squiu umigc is (loliiicd vs 


u(x,?/) 


1 if |a:| < 0 5 , |y| <05 
0 otherwise 


(6 1) 


1 ho piojcctioiib of the square hvvc a trapezoidal shape as discussed in appendix A and 


can bo calculated in a straightforward manner The projections are sampled so as to get 
200 samples pc r projceLion, ic M = 200 and no of projection angles N = 128 The square 
mngc IS icconsliurtecl on a giid of 100 x 100 This gives a oversampling ratio p = 2 
llu iccoiisliuction IS cnuiccl out without idclmg any noise to the projoctions DilTcront 
convolve IS u< use d loi compulsion purpose Simil vtly, using 400 samples per projections, 


which gives ovcisainplmg latio p = 4, reconstruction is earned out The reconstructed 
imvgc IS noinnh/(d so is to get a value of 1 0 at the center pixel 

Wc (O.npuc the pcilo.m.iKC ofchfrccnl convolves let u(A , /) dc note the ic c oust . uc t c cl 
imige attci normalization We define tlie Signal to Noise ratio SNR as 


5iVi? = lOln 




(6 2 ) 
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I he ic Hulls lie pit sc 


iil( d 111 I d)l( 0 1 



MSL 



1 ijjUK () 1 MSI Vs M \silh (onsl iiiL N 

I Ills I ibl( shows iliil foi Ilni filld lli< SNR inijirovc s wilh mcic \sni{, vRik of (> while 
loi oIIki filleis the pcifoiininec dc f( iioiiak s incl llir SNR fills Iliis usulls liowevei 
docs tiol itnfch willi flu (iioi (uivt j,ivcn iii [25] Ihc cutve is uiuodiicecl Ik u in li^ 
() I I hus I ho Mean Squ u( I not is shown to siluratc to i nmiiniuni vihu ilfi i A/ 
c'{((((ls in opfiimim vduc {jivon by M = 2 * A^/tt Wc foci tint Ihis disciciniuy could 
be due lo me leased qu ml i/ ilioii ,eom|)uhli<)ii uid loimelmg noise winch mcie ise s with 
me i( ise el h md width Alsei with me le ise d h viidwidlh the ic ( onsli m Ic el squ iic h is liij^lic i 
oveishool md uiidcishool doiif, Iho odgos IIowovoi for llaii lillci the ic eoiisl i in I loii 
lilleiis i futic I lou ol ove IS uiipliii^ 1 U lo p Aiid it L decs into iccouiil the meic ise d ptop 1 1 uni 

b mdwidi h 

Nik Hut the Run 1 ik fillci Ins nmp liko lesponsc with simp cutoff 11ns lesulls in 
simp ie uisitioei e^f the squ uo edges with high overshoot and ringing ihc oveishool ind 
imgmg e in be le eluced by h wmg isinooth roll off as in case of Gene i ili/e el II uiiimiig Idle i 
llowcvei Hub is at Hio cost of wider trinsition width llie p ii imotoi ev of the C II fillci 

blows 1 (I ideotr be twee 11 oveishool -uiel liuisitioii wiclfli rv = 1 e one sponels lo R im I ik 
hlkl loi compulsion pm pose we hive seleeleel leeiilei Vilue ol I) 5 1 he e leiss se e I ion 
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i iblo 6 1 Performance comparison (No noise added) 
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Figuie 6 11 Phantom reconstruction using Generalized Hamming filter /? = 4 



1 iguic () 12 Phantom reconstruction using Ilaai filter p - i 
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6 2 Reconstruction from Noisy projections 


Now wc con&idci Uic rtcoustiuction from noisy projections using wavelet preprocessing 
I he (ill Cl used loi wivclct dccoiuposition is Daubcclucs 8 (ilic bimulatioii was also car 
lied out with filtcis of length 4, 10, and 20 However there was not significant change in 
pcifoiiintico ) 

Lhc viluc ol bNH IS evduated foi diflcicut couvolvcis at diiferciit ovei sampling ratio p and 
diffeient noise poweis Two cases are considered separately In first case white gaussian 
noise of /jcio me XU is added to the piojcctions No wavelet pre processing is carried and 
the SNR IS mcasuicd foi diffeient noise powers These results are suimnarized m table 6 2 
and 6 4 In the second case, pre processing of noisy projections, as discussed in chp5 is 
can led out befoio convolution and backprojection The results are summarized in table 
6 3 ihc same pioceduie is repeated for p = 2 and the results are presented in table 6 3 
uid 6 5 In dl I lie sc eases, tlic iioimahzatioii factor used is same as m case of noiseless 
lecoiistiuction with respective filters and oversampling ratio 
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iable 6 2 Perioiinance comparision at various noise powers No pre processing p — 4 
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1 ible 6 3 Peiloimxnce compaiisioii at various noise poweis, with wavelet pro processing 

p = 4 
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Figure 6 17 Noisy phantom reconstruction using wavelet pre processing and Generalized 
Ilamnuiig filtei p — A 



1 iguie 6 18 Noisy phantom reconstruction using wavelet pre processing and Haar filter 



Now wc conside: the peiforraance with constant MJN ratio le the no of projection 
ingles aie incic xsed in piopoition to the sampling rate in projection domain The SNR 
IS mcasuicd loi the squue image as well as for the Shepp Logan phantom For the square 
image ,thc II i u filtei shows betlei peiformance For Shepp Logan phantom, the Llaar filter 
again gives betici icsuHs in compaiision with the Ram Lak filter However the performance 
ol Ram I ak impioves with increasing iV, while Ilaai filter has shown a small decrease m 
the SNR value These lesults aie sumanzed in table 6 7 &: table 6 8 
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1 ibk 6 7 PciloimaiKc compaiision with constant NjM ratio (Using square image) 
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libk 6 8 Ikiloimaucc compulsion with constant NjM ratio (Using Shepp Logan 
phantom) 


6 3 Conclusion 

I„ tius wolk W( have lUcmplM lo study the relative peitormauco of different filters used 
in imige rcionstulction the mam limitation of Ham filter rs that it is not defined for 
odd vJms ol oYcisinipling ratio The Ilaar filter, however, has given improved results 
for certain .lass of imiges such as square images We feel that tor .mages belonging 
to a moic geneiil class, additional simulation work should he earned out m order to 
measure 1 he iclai ive performance of Ilaar filter Also the Haar filter wavelet pre processing 
combrnation gives better results These conclusions are based on the results of compu er 
smuUtions and not based on any kind of analytical treatment 
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Appendix A 


ilic Wq € [0,7r/4] has a geneial trapezoidal shape as shown in hg A 1 where 

using simple tnngnometiy we can show that 


Ofl = — (cos0 + sm0) (A* 1) 

Ct 

bs = ^{cosd — smO) (^ 2 ) 


It IS shown 111 [14] (hat 



ir;(kl) = " 


Jl_ 


(A3) 
(A 4) 


Hence fiom cq 3 20 wc get 

k(tJ) = -^*wot,(t,e) (A 6) 

roi the sake of simplicity wc will drop the term 0 m the discussion to follow and we also 
assume that ko — 1 Hence 



Figure A 1 Geometry of uio a(-s, 0) 
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In 


(A 8) 


[lnlt-rl]j (A 9) 


(A 10) 
(All) 
(A 12) 

(rt~6) \l^~b^\ 

Aftei subslitulin^ Uie Vdlues of ag and bg as given by eq A 1 , A 2 and myltiplymg by hg 
we get the final result as 


1 \t^ - -y(l 4-sm2^)| 

” 7rsin20 |t2 - M(l - sin 20)1 

The above discussion is valid for 0 0 Ihe expression for A:(t,0) can 

limiting case of the above equation and the reult is 


(A 13) 

be evaluated as a 




(A 14) 
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